Ce script décrit plusieurs fonctions de manipulation de la dimension spatiale des données DVF. Il propose aussi plusieurs formes de représentation cartographique de ces données à une échelle de l’EPCI (La Rochelle agglomération).
library(tidyverse) # Manipulation de données
library(sf) # Manipulation de données spatiales
library(mapsf) # Cartographie thématique
library(cartography) # Cartographie thématique
library(linemap) # cartes en pics
library(spdep) # pour calculer l'auto-corrélation spatiale
library(geoR) # pour calculer le semi-variogramme empirique
library(spatstat) # pour produire des surfaces lissées
library(terra) # pour le traitement de données matricielles (raster)
library(cluster) # pour la CAH
library(happign) #pour récupérer les référentiels de l'IGN
library(maptiles) #pour ajouter des fonds aux représentattions carto
library(stars) #pour convertir des objets ppp
library(raster) #pour découper un raster par l'emprise d'un vecteur
DVF <- read.csv("data/dvf_rochelle.csv", encoding = "UTF-8")
RecapLR <- DVF %>% group_by(type) %>%
summarise(tot = n(), prixmed = round(median(prix)), prixmoy = round(mean(prix)), surfmed = round(median(surface)), surfmoy = round(mean(surface)), prixm2med = round(median(prix_m2)), prixm2moy = round(mean(prix_m2)))
RecapLR <- RecapLR %>% mutate(part = (tot/sum(tot)*100))
print(RecapLR)
## # A tibble: 2 × 9
## type tot prixmed prixmoy surfmed surfmoy prixm2med prixm2moy part
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Appartement 9377 145000 174955 45 50 3562 3625 36.5
## 2 Maison 16279 247000 279364 95 102 2652 2806 63.5
cog <- happign::get_wfs(
x = NULL,
apikey = "administratif",
layer = "ADMINEXPRESS-COG-CARTO.LATEST:commune",
ecql_filter = "insee_dep = '17'"
)
## Features downloaded : 463
epci_lr <- cog %>%
filter(siren_epci == "241700434")
# Récupération des Iris par la liste des communes de l'EPCI
liste_com_epci <- list(unique(epci_lr$insee_com))
liste_com_epci <- stringi::stri_join_list(liste_com_epci, sep="','")
IRISLR <- happign::get_wfs(
x = NULL,
apikey = "cartovecto",
layer = "STATISTICALUNITS.IRIS:contours_iris",
ecql_filter = paste0('insee_com in (\'', liste_com_epci, '\')')
)
## Features downloaded : 67
# Récupération des communes de l'EPCI et des communes attenantes
cog$inter <- st_intersects(cog,st_union(epci_lr), sparse=F)
com_LR_CA <- cog %>% filter(inter == T)
## Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
## ℹ Please use one dimensional logical vectors instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
liste_com_lrca <- list(unique(com_LR_CA$insee_com))
liste_com_lrca <- stringi::stri_join_list(liste_com_lrca, sep="','")
IRISLR_CA <- happign::get_wfs(
x = NULL,
apikey = "cartovecto",
layer = "STATISTICALUNITS.IRIS:contours_iris",
ecql_filter = paste0('insee_com in (\'', liste_com_lrca, '\')')
)
## Features downloaded : 80
cog <- cog %>% st_transform(2154)
com_LR_CA <- com_LR_CA %>% st_transform(2154)
IRISLR <- IRISLR %>% st_transform(2154)
IRISLR_CA <- IRISLR_CA %>% st_transform(2154)
# Calcul du prix moyen au m² par IRIS dans La Rochelle Agglomération
DVF <- st_as_sf(DVF, coords = c('longitude', 'latitude'), crs = 4326)
DVF <- st_transform(DVF, 2154)
IRISDVF <- aggregate(dplyr::select(DVF, prix_m2), by = IRISLR, FUN = mean)
IRISDVF <- IRISDVF %>% drop_na()
# Sélection des communes dont le nom sera affiché sur la carte (top 6 selon le nb d'habitants)
liste_noms_LR <- c("La Rochelle", "Aytré", "Puilboreau", "Châtelaillon-Plage", "Lagord", "Périgny")
Selection_Communes_LR <- cog[which(cog$nom %in% liste_noms_LR),]
# Affichage de la carte
par(oma = c(1,0,0,0))
# Génération du fond de carte
fd1 <- get_tiles(IRISDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
plot.new()
mf_raster(fd1)
# Génération de la carte choroplète
mf_map(
type = "choro",
x = IRISDVF,
var = "prix_m2",
method = "quantiles",
nbreaks = 5,
pal = "YlOrRd",
border = "black",
leg_title = "Prix moyen \nau m² (Euros)",
leg_pos = c(399000,6564450),
add = TRUE)
## The following argument is not relevant when using type = 'choro': method.
# Ajouter les étiquettes des communes selectionnées plus haut, et les autres éléments d'habillage
mf_label(
Selection_Communes_LR,
var = "nom",
halo = TRUE,
bg = "white",
r = 0.1,
cex = 0.8,
overlap = FALSE)
mf_title("Prix moyen au m² de l'immobilier par Iris sur La Rochelle agglomération et communes attenantes (2014-2022)", cex = 0.8)
mf_credits("IGN et DGFip")
mf_scale(size = 5)
mf_arrow("topleft")
# Calcul du nombre de transactions par IRIS
IRISDVF$cpt_vente <- lengths(st_intersects(IRISDVF,DVF))
IRISDVF <- IRISDVF %>% drop_na()
# Affichage de la carte
par(oma = c(1,0,0,0))
# Génération du fond de carte
fd1 <- get_tiles(IRISDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
plot.new()
mf_raster(fd1)
plot(st_geometry(IRISDVF), add=TRUE)
# Génération de la carte en symboles proportionnels
mf_map(
type = "prop",
x = IRISDVF,
var = "cpt_vente",
inches = 0.15,
col = rgb(0,0.2,1,alpha = 0.5),
border = rgb(0,0,0,alpha = 0.5),
leg_title = "Nombre de transactions \npar iris",
add = TRUE)
# Ajouter les étiquettes des communes selectionnées plus haut, et les autres éléments d'habillage
mf_label(
Selection_Communes_LR,
var = "nom",
halo = TRUE,
bg = "white",
r = 0.1,
cex = 0.8,
overlap = FALSE)
mf_title("Nombre de mutations par Iris (2014-2022) à La Rochelle agglomération et communes attenantes", cex = 0.8)
mf_credits("IGN et DGFip")
mf_scale(size = 5)
mf_arrow("topleft")
# Récupération via une liste des communes de l'EPCI
liste_com <- unique(IRISLR$insee_com)
liste_com <- gsub("^(\\d{2})(.*)", "\\2", liste_com, perl=T)
st <- list(liste_com)
str_com <- stringi::stri_join_list(st, sep="','")
sec_cad <- happign::get_wfs(
x = NULL,
apikey = "parcellaire",
layer = "BDPARCELLAIRE-VECTEUR_WLD_BDD_WGS84G:divcad",
ecql_filter = paste0('code_com in (\'', str_com, "') and code_dep = '17'")
)
## Features downloaded : 564
# Même opération en rajoutant les communes attenantes à l'EPCI
liste_com_ca <- unique(IRISLR_CA$insee_com)
liste_com_ca <- gsub("^(\\d{2})(.*)", "\\2", liste_com_ca, perl=T)
st <- list(liste_com_ca)
str_com <- stringi::stri_join_list(st, sep="','")
sec_cad_ca <- happign::get_wfs(
x = NULL,
apikey = "parcellaire",
layer = "BDPARCELLAIRE-VECTEUR_WLD_BDD_WGS84G:divcad",
ecql_filter = paste0('code_com in (\'', str_com, "') and code_dep = '17'")
)
## Features downloaded : 848
plot(sec_cad["code_com"])
sec_cad <- st_transform(sec_cad,2154)
sec_cad_ca <- st_transform(sec_cad_ca,2154)
# Calcul du prix moyen au m² par section
sec_cadDVF <- sec_cad %>%
st_join(DVF) %>%
group_by(id) %>%
summarise(prix_m2 = round(mean(prix_m2))) %>%
na.omit()
# Affichage de la carte
par(mar = c(1, 0, 0, 0))
# Génération du fond de carte
fd1 <- get_tiles(IRISDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
plot.new()
mf_raster(fd1)
plot(st_geometry(sec_cad), add=TRUE)
# Génération de la carte choroplète
mf_map(
type = "choro",
x = sec_cadDVF,
var = "prix_m2",
method = "quantiles",
nbreaks = 6,
pal = "YlOrRd",
border = rgb(0,0,0,alpha=0.5),
leg_title = "Prix moyen \nau m² (Euros)",
add = TRUE)
## The following argument is not relevant when using type = 'choro': method.
# Ajout des étiquettes et autres éléments d'habillage
mf_label(
Selection_Communes_LR,
var = "nom",
halo = TRUE,
bg = "white",
r = 0.1,
cex = 0.8,
overlap = FALSE)
mf_title("Prix moyen au m² de l'immobilier par section cadastrale à La Rochelle agglo (2014-2022)", cex = 0.8)
mf_credits("IGN et DGFip")
mf_scale(size = 5)
mf_arrow("topleft")
# Calcul du nombre de mutations par section
sec_cadDVF <- sec_cadDVF %>%
mutate(Nbmutations = lengths(st_intersects(sec_cadDVF, DVF)))
# Affichage de la carte
par(mar = c(1, 0, 0, 0))
# Génération du fond de carte
fd1 <- get_tiles(IRISDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
mf_raster(fd1)
plot(st_geometry(sec_cad), border= rgb(0,0,0,alpha=0.3), add=T)
# Génération de la carte en symboles proportionnels
propSymbolsLayer(x = sec_cadDVF,
var = "Nbmutations",
col = rgb(0,0.2,1,alpha = 0.5),
border = "white",
lwd = 0.01,
inches = 0.1,
add = TRUE,
legend.pos = "n")
# Ajout des étiquettes et autres éléments d'habillage
layoutLayer(title = "Nombre de mutations par section cadastrale à La Rochelle Agglomération (2014-2022)", source = " IGN et DGFip", north = TRUE, horiz = FALSE, tabtitle = TRUE, frame = FALSE, col = "black", coltitle = "white")
legendCirclesSymbols(
pos = "bottomleft",
title.txt = "Nombre de mutations DVF",
title.cex = 0.8,
cex = 1,
border = "white",
lwd = 0.01,
values.cex = 0.6,
var = c(min(sec_cadDVF$Nbmutations), max(sec_cadDVF$Nbmutations)),
inches = 0.08,
col = rgb(0,0.2,1,alpha = 0.5),
frame = FALSE,
values.rnd = 0,
style = "e"
)
# Ajouter les étiquettes des communes selectionnées plus haut
labelLayer(Selection_Communes_LR, txt = "nom", halo = TRUE, bg = "white", r = 0.1, cex = 0.8, pos = 3, font = 3, offset = 0, overlap = TRUE)
# pour calculer les pseudo-centroides des Iris
IRISLR.Centroids <- st_point_on_surface(IRISLR)
# Identification du plus proche voisin de chaque Iris
listPPV_iris_cen <- knearneigh(IRISLR.Centroids, k = 1)
# Conversion de l'objet knn en objet nb
PPV_iris <- knn2nb(listPPV_iris_cen, row.names = IRISLR.Centroids$code_iris)
distPPV_iris <- nbdists(PPV_iris, IRISLR.Centroids) # distance entre plus proches voisins
print(as.data.frame(t(as.matrix(summary(unlist(distPPV_iris)))))) # résumé statistique
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 331.1623 543.4473 1081.793 1339.615 1882.393 4854.631
# des distances entre IRIS
hist(unlist(distPPV_iris), breaks = 20,
main = "Distance au plus proche voisin",
col = "black", border = "white", xlab = "Distance", ylab = "Fréquence")
# pour convertir les Iris en objet nb
IRIS.nb <- poly2nb(pl = IRISDVF,
snap = 50,
queen = TRUE) # voisinage au sens large sommet ou frontière communes
#calcul du test de Moran
moran.test(IRISDVF$prix_m2, listw = nb2listw(IRIS.nb))
##
## Moran I test under randomisation
##
## data: IRISDVF$prix_m2
## weights: nb2listw(IRIS.nb)
##
## Moran I statistic standard deviate = 8.1923, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.620170543 -0.015151515 0.006014213
Relation supérieur (Moran I) à 0 doc autocorrélation positive. Z-score = 8.1923, en dehors de l’intervalle [-1,96 ; 1,96] pour \(\alpha\) = 0.95~%. On lisse quand on a de l’autocorrélation spatiale, sinon c’est déconseillé. Ici, l’indice de Moran nous autorise à représenter les données par un lissage spatial
Nous décidons de réaliser le lissage en prenant en compte les communes attenantes afin de limiter l’effet de bord.
# pour définir le contour de La Rochelle agglomération (et communes attenantes) comme emprise pour le lissage
EmpriseLR <- as.owin(st_geometry(com_LR_CA))
IRISDVF_CA <- aggregate(dplyr::select(DVF, prix_m2), by = IRISLR_CA, FUN = mean)
IRISDVF_CA <- IRISDVF_CA %>% drop_na()
# pour récupérer les coordonnées du centroïde des IRIS
pts <- st_coordinates(st_point_on_surface(st_geometry(IRISLR_CA)))
# pour créer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (prix au m²)
IRISDVF_CA.ppp <- ppp(pts[,1], pts[,2], window = EmpriseLR, marks = IRISDVF_CA$prix_m2)
# pour calculer la surface lissée (rayon lissage : 1 km et resolution spatiale de l'image : 1 ha)
cartelissee <- Smooth(IRISDVF_CA.ppp, sigma = 1000, weights = IRISDVF_CA.ppp$marks, eps = 100)
# Conversion de la surface lissée au format raster et découpage selon l'emprise
# de l'EPCI
# Conversion de la surface lissée au format raster
cartelissee.raster <- terra::rast(cartelissee)
crs(cartelissee.raster) <- "epsg:2154"
# Convertir epci_lr en spatVector
mask <- vect(st_transform(epci_lr,2154))
masked <- terra::mask(cartelissee.raster, mask)
# Paramétrage des marges de la fenêtre pour maximiser l'emprise de la carte
par(mar = c(1, 0, 0, 0))
# Pour afficher la surface lissée et définir l'habillage
# Calcul des seuils
bks <- unique(getBreaks(values(masked), method = "q6"))
# Création d'une palette de couleurs (double gradation harmonique)
cols <- c("#2a9c4e", "#77c35c", "#c4e687", "#ffffc0", "#fec981", "#dc292c")
fd1 <- get_tiles(IRISLR_CA,"CartoDB.PositronNoLabels", 12, crop=T)
fd1 <- project(fd1, "epsg:2154")
mf_raster(fd1)
# Affichage de la carte lissée et du contour des Iris
plot(masked, breaks = bks, col = cols, add = T, legend = F)
plot(epci_lr$geometry, border = rgb(1,1,1, alpha = 0.5), lwd = 0.5, add = T)
plot(IRISDVF$geometry, border = "grey60", lwd = 0.05, lty = 3, add = T)
legendChoro(
pos = "bottomleft",
title.txt = "Prix moyen \nau m² (Euros)",
breaks = bks,
nodata = FALSE,
values.rnd = -1,
border = "white",
col = cols
)
# Habillage de la carte
layoutLayer(title = "Prix moyen au m² de l'immobilier à La Rochelle Agglo (2014-2022)",
author = " Sources : IGN et DGFip - Rayon : 1000 m,\n résolution : 100m - A partir des Iris",
scale = 5, frame = TRUE, col = "black", coltitle = "white",
north = TRUE, tabtitle = TRUE, horiz = FALSE)
# Ajouter les étiquettes des communes sélectionnées plus haut
labelLayer(Selection_Communes_LR, txt = "nom", halo = TRUE, bg = "white", r = 0.05, cex = 0.6, pos = 3, font = 3, offset = 0, overlap=F)
# Récupération des valeurs de prix_m2 par section cadastrale
sec_cad_caDVF <- sec_cad_ca %>%
st_join(DVF) %>%
group_by(id) %>%
summarise(prix_m2 = round(mean(prix_m2)))
sec_cad_caDVF <- sec_cad_caDVF %>% na.omit()
# pour calculer les pseudo-centroides des sections cadastrales (indispensable pour le calcul sur les plus proches voisins)
sec_cad.Centroids <- st_point_on_surface(sec_cad_caDVF)
# Calcul sur les plus proches voisins
listPPV <- knearneigh(sec_cad.Centroids, k = 1) # pour connaître le plus proche voisin de chaque section cadastrale
PPV <- knn2nb(listPPV, row.names = sec_cad.Centroids$id) # pour convertir l'objet knn en objet nb
distPPV <- nbdists(PPV, sec_cad.Centroids) # pour connaître la distance entre plus proches voisins
print(as.data.frame(t(as.matrix(summary(unlist(distPPV))))))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 126.1835 405.637 514.648 581.8634 712.9702 3162.455
hist(unlist(distPPV), breaks = 20,
main = "Distance au plus proche voisin",
col = "black", border = "white", xlab = "Distance", ylab = "Fréquence")
# pour convertir les sections en objet nb
sec_cad_ca.nb <- poly2nb(pl = sec_cad_caDVF,
snap = 980,
queen = TRUE) #Nous augmentons le snap à 980, car l'indice
#de Moran ne peut pas être calculé avec une valeur inférieure
#calcul du test de Moran
moran.test(sec_cad_caDVF$prix_m2, listw = nb2listw(sec_cad_ca.nb))
##
## Moran I test under randomisation
##
## data: sec_cad_caDVF$prix_m2
## weights: nb2listw(sec_cad_ca.nb)
##
## Moran I statistic standard deviate = 35.789, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.6004203228 -0.0016977929 0.0002830445
# pour récupérer les coordonnées du centroïde des sections
pts <- st_coordinates(st_point_on_surface(st_geometry(sec_cad_caDVF)))
# pour créer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (prix au m²)
sec_cad.ppp <- ppp(pts[,1], pts[,2], window = EmpriseLR, marks = sec_cad_caDVF$prix_m2)
## Warning: 4 points were rejected as lying outside the specified window
# pour calculer la surface lissée (rayon lissage : 1 km et resolution spatiale de l'image : 1 ha)
cartelissee <- Smooth(sec_cad.ppp, sigma = 1000, weights = sec_cad.ppp$marks, eps = 100)
# Conversion de la surface lissée au format raster
cartelissee.raster <- terra::rast(cartelissee)
crs(cartelissee.raster) <- "epsg:2154"
# Convertir epci_lr en spatVector
mask <- vect(st_transform(epci_lr,2154))
masked <- terra::mask(cartelissee.raster, mask)
# Paramétrage des marges de la fenêtre pour maximiser l'emprise de la carte
par(mar = c(1, 0, 0, 0))
# pour afficher la surface lissée et définir l'habillage
# Calcul des seuils
bks <- unique(getBreaks(values(masked), method = "q6"))
# Création d'une palette de couleurs (double gradation harmonique)
cols <- c("#2a9c4e", "#77c35c", "#c4e687", "#ffffc0", "#fec981", "#dc292c")
fd1 <- get_tiles(sec_cad_caDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
mf_raster(fd1)
# Affichage de la carte lissée et du contour des Iris
plot(masked, breaks = bks, col = cols, add = T, legend = F)
plot(sec_cad_ca$geometry, border = rgb(0,0,0,alpha = 0.5), lwd = 0.03, lty = 3, add = T)
plot(com_LR_CA$geometry, border = "white", lwd = 0.5, add = T)
legendChoro(
pos = "bottomleft",
title.txt = "Prix moyen \nau m² (Euros)",
breaks = bks,
nodata = FALSE,
values.rnd = -1,
border = "white",
col = cols
)
# Habillage de la carte
layoutLayer(title = "Prix moyen au m² de l'immobilier à La Rochelle Agglo (2014-2022)",
author = " Sources : IGN et DGFip - Rayon : 1000 m,\n résolution : 100m - A partir des Iris",
scale = 5, frame = TRUE, col = "black", coltitle = "white",
north = TRUE, tabtitle = TRUE, horiz = FALSE)
# Ajouter les étiquettes des communes sélectionnées plus haut
labelLayer(Selection_Communes_LR, txt = "nom", halo = TRUE, bg = "white", r = 0.05, cex = 0.6, pos = 3, font = 3, offset = 0, overlap=F)
# Calcul sur les plus proches voisins
listPPV <- knearneigh(DVF, k = 1) #pour connaître le plus proche voisin de chaque mutation
PPV <- knn2nb(listPPV, row.names = DVF$id_vente) # pour convertir l'objet knn en objet nb
distPPV <- nbdists(PPV, DVF) # pour connaître la distance entre plus proches voisins
print(as.data.frame(t(as.matrix(summary(unlist(distPPV))))))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 0 0 6.784495 17.69769 26.6294 2459.984
hist(unlist(distPPV), nclass = 500,
main = "Distance au plus proche voisin",
col = "black", border = "white", xlab = "Distance", ylab = "Fréquence", xlim = c(0,150))
#calcul du test de Moran
moran.test(DVF$prix_m2, listw = nb2listw(PPV))
##
## Moran I test under randomisation
##
## data: DVF$prix_m2
## weights: nb2listw(PPV)
##
## Moran I statistic standard deviate = 61.473, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 4.690504e-01 -3.897876e-05 5.823001e-05
# Extraction des coordonnées spatiales
coordinates <- st_coordinates(DVF)
# pour créer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (prix moyen au m²)
DVF.ppp <- ppp(x = coordinates[, 1], y = coordinates[, 2], window = EmpriseLR, marks = DVF$prix_m2)
# pour calculer la surface lissée (rayon lissage : 0,5 km et résolution spatiale de l'image : 1 ha) --> calcul long
cartelissee <- Smooth(DVF.ppp, sigma = 500, weights = DVF.ppp$marks, eps = 100)
# Conversion de la surface lissée au format raster
cartelissee.raster <- terra::rast(cartelissee)
crs(cartelissee.raster) <- st_crs(DVF)$srid
# découpage de la surface lissée sur l'emprise Rennes Métropole
cartelissee.raster <- mask(cartelissee.raster, st_transform(epci_lr,2154))
# Paramétrage des marges de la fenêtre pour maximiser l'emprise de la carte
par(mar = c(1, 0, 0, 0))
# pour afficher la surface lissée et définir l'habillage de la carte
# Calcul des seuils
bks <- unique(getBreaks(values(cartelissee.raster), method = "q6"))
# Création d'une palette de couleurs (double gradation harmonique)
cols <- c("#2a9c4e", "#77c35c", "#c4e687", "#ffffc0", "#fec981", "#dc292c")
fd1 <- get_tiles(DVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
mf_raster(fd1)
# Affichage de la carte lissée et du contour des Iris
plot(cartelissee.raster, breaks = bks, col = cols, add = T, legend = F)
plot(com_LR_CA$geometry, border = "white", lwd = 0.5, add = T)
legendChoro(
pos = "bottomleft",
title.txt = "Prix moyen \nau m² (Euros)",
breaks = bks,
nodata = FALSE,
values.rnd = -1,
border = "white",
col = cols
)
# Habillage de la carte
layoutLayer(title = "Prix moyen au m² de l'immobilier à La Rochelle Agglo (2014-2022)",
author = " Sources : IGN et DGFip - Rayon : 500 m,\n résolution : 100m - A partir des Iris",
scale = 5, frame = TRUE, col = "black", coltitle = "white",
north = TRUE, tabtitle = TRUE, horiz = FALSE)
# Ajouter les étiquettes des communes sélectionnées plus haut
labelLayer(Selection_Communes_LR, txt = "nom", halo = TRUE, bg = "white", r = 0.05, cex = 0.6, pos = 3, font = 3, offset = 0, overlap=F)
Une typologie est établie à partir d’indicateurs immobiliers et une cartographie des sous-marchés associés est réalisée.
MutationsIRIS <- DVF %>%
st_join(IRISLR) %>%
drop_na()
MutationsIRIS <- as.data.frame(MutationsIRIS)
MutationsIRIS <- MutationsIRIS %>%
mutate(prixm2 = prix / surface)
IRISDVFClassif1 <- MutationsIRIS %>%
group_by(code_iris) %>%
summarise(Nbtransactions = n(),
code_iris = first(code_iris),
Prixmoyen = round(mean(prix)),
Prixm2moyen = round(mean(prixm2)),
Surfacemoyenne = round(mean(surface)),
PropMaison = length(type[type == "Maison"]) / Nbtransactions * 100,
PropAppart = length(type[type == "Appartement"]) / Nbtransactions * 100)
IRISDVFClassif <- data.frame(IRISDVFClassif1[, c("Nbtransactions", "Prixmoyen", "Prixm2moyen", "Surfacemoyenne", "PropMaison", "PropAppart")])
Les variables sont centrées réduites afin de standardiser les données.
IRISDVFClassifScale <- scale(IRISDVFClassif)
CAHIRIS <- agnes(
IRISDVFClassifScale,
metric = "euclidean",
method = "ward"
)
sortedHeight <- sort(CAHIRIS$height, decreasing = TRUE)
relHeight <- sortedHeight/ sum(sortedHeight) * 100
cumHeight <- cumsum(relHeight)
library(Factoshiny)
## Le chargement a nécessité le package : FactoMineR
## Le chargement a nécessité le package : shiny
## Le chargement a nécessité le package : FactoInvestigate
library(FactoMineR)
pca <- PCA(IRISDVFClassifScale, ncp = Inf, scale.unit = FALSE, graph = FALSE)
cah <- HCPC(pca, graph = FALSE)
plot(cah, choice = "bar")
# Gain d'inertie
barplot(
relHeight[1:30],
names.arg = seq(1, 30, 1),
col = "black", border = "white",
xlab = "Noeuds", ylab = "Part de l'inertie totale (%)"
)
barplot(
cumHeight[1:30],
names.arg = seq(1, 30, 1),
col = "black", border = "white",
xlab = "Nombre de classes", ylab = "Part de l'inertie totale (%)"
)
dendroCSP <- as.dendrogram(CAHIRIS)
plot(dendroCSP, leaflab = "none")
clusIRIS <- cutree(CAHIRIS, k = 3)
IRISCluster <- IRISDVFClassif1
IRISCluster$CLUSIMMO <- factor(
clusIRIS,
levels = 1:3,
labels = paste("Classe", 1:3)
)
RecapCAHIRIS <- IRISCluster %>%
group_by(CLUSIMMO) %>%
summarise(
NB= n(),
NbTransac = round(mean(Nbtransactions)),
Prixmoyen = round(mean(Prixmoyen)),
Prixm2 = round(mean(Prixm2moyen)),
Surface = round(mean(Surfacemoyenne)),
PropMaison = mean(PropMaison),
PropAppart = mean(PropAppart)
)
print(RecapCAHIRIS)
## # A tibble: 3 × 8
## CLUSIMMO NB NbTransac Prixmoyen Prixm2 Surface PropMaison PropAppart
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Classe 1 37 293 271011 2722 102 90.1 9.92
## 2 Classe 2 18 182 197212 2611 77 66.9 33.1
## 3 Classe 3 12 713 217482 3775 60 23.7 76.3
SyntheseCAHIRIS <- RecapCAHIRIS %>% mutate(
nbtransacmoy = mean(IRISDVFClassif$Nbtransactions),
surfacemoy = mean(IRISDVFClassif$Surfacemoyenne),
prixmoy = mean(IRISDVFClassif$Prixmoyen),
prixm2moyen = mean(IRISDVFClassif$Prixm2moyen),
propmaisonmoyen = mean(IRISDVFClassif$PropMaison),
propappartmoyen = mean(IRISDVFClassif$PropAppart),
NbMutations = (NbTransac - nbtransacmoy) / nbtransacmoy * 100,
Prix = (Prixmoyen - prixmoy) / prixmoy * 100,
Prixm2 = (Prixm2 - prixm2moyen) / prixm2moyen * 100,
Surface = (Surface - surfacemoy) / surfacemoy * 100,
PropMaison = (PropMaison - propmaisonmoyen) / propmaisonmoyen * 100,
PropAppart = (PropAppart - propappartmoyen) / propappartmoyen * 100
)
# pour créer un nouveau tableau avec une sélection de colonnes à faire apparaître sur le graphique
SyntheseCAHIRIS <- data.frame(
SyntheseCAHIRIS[,
c("CLUSIMMO", "NbMutations", "Surface", "Prix", "Prixm2", "PropMaison", "PropAppart")
]
)
gather <- SyntheseCAHIRIS %>%
gather(key = variable, value = "value", NbMutations:PropAppart)
ggplot(gather, aes(x = variable, y = value, fill = CLUSIMMO)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_manual(values = c("#2A9D8F","#264653","#fcbf49","#e45c3a","#0074D9")) +
ylab("Variation par rapport à la moyenne métropolitaine (%)") +
theme_bw() +
theme(legend.position = "none") +
facet_wrap(~CLUSIMMO, ncol = 1)
# pour joindre le résultat de la typologie dans la couche des IRIS
IRISDVFCAH <- left_join(IRISLR, IRISCluster, by = "code_iris")
par(mar = c(0, 0, 1.2, 2))
# Affichage des communes de RM et alentours en arrière-plan et centrage de la carte sur Rennes Métropole
plot(
st_geometry(IRISLR),
lwd = 0.1,
col = "grey95",
border = "white",
bg = "#B5D0D0"
# xlim = c(st_bbox(IRISRM)[1], st_bbox(IRISRM)[3]), ylim = c(st_bbox(IRISRM)[2], st_bbox(IRISRM)[4])
)
# Affichage de la typologie et du contour des IRIS
unique(IRISDVFCAH$CLUSIMMO)
## [1] Classe 1 Classe 2 Classe 3
## Levels: Classe 1 Classe 2 Classe 3
IRISDVFCAH
## Simple feature collection with 67 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 372869.8 ymin: 6552942 xmax: 398644.9 ymax: 6581807
## Projected CRS: RGF93 v1 / Lambert-93
## # A tibble: 67 × 17
## id insee_com nom_com iris code_iris nom_iris typ_iris layer path
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 IRIS____0000… 17291 Puilbo… 0000 172910000 Puilbor… Z CONT… R:/c…
## 2 IRIS____0000… 17300 La Roc… 0505 173000505 Mireuil… H CONT… R:/c…
## 3 IRIS____0000… 17300 La Roc… 0202 173000202 La Gene… H CONT… R:/c…
## 4 IRIS____0000… 17274 Périgny 0103 172740103 Périgny… H CONT… R:/c…
## 5 IRIS____0000… 17420 Salles… 0000 174200000 Salles-… Z CONT… R:/c…
## 6 IRIS____0000… 17300 La Roc… 0402 173000402 La Pall… H CONT… R:/c…
## 7 IRIS____0000… 17153 Esnand… 0000 171530000 Esnandes Z CONT… R:/c…
## 8 IRIS____0000… 17300 La Roc… 1001 173001001 Les Min… H CONT… R:/c…
## 9 IRIS____0000… 17200 Lagord 0103 172000103 Ouest H CONT… R:/c…
## 10 IRIS____0000… 17300 La Roc… 0903 173000903 Villene… H CONT… R:/c…
## # ℹ 57 more rows
## # ℹ 8 more variables: geometry <MULTIPOLYGON [m]>, Nbtransactions <int>,
## # Prixmoyen <dbl>, Prixm2moyen <dbl>, Surfacemoyenne <dbl>, PropMaison <dbl>,
## # PropAppart <dbl>, CLUSIMMO <fct>
typoLayer(
x = IRISDVFCAH,
var="CLUSIMMO",
col = c("#2A9D8F","#264653","#fcbf49","#e45c3a","#0074D9"),
lwd = 0.05,
border = "grey70",
legend.pos = "bottomleft",
legend.title.txt = "Sous-marchés \nimmobiliers",
legend.nodata = "Aucune mutation",
add = TRUE)
layoutLayer(title = "Sous-marchés immobiliers dans Rennes Métropole à l'échelon des IRIS (2014-2019)",
author = " Sources : IGN et DGFip - Typologie obtenue par CAH",
scale = 5, frame = TRUE, col = "black", coltitle = "white",
north = TRUE, tabtitle = TRUE, horiz = FALSE)
# Ajouter les étiquettes des communes sélectionnées plus haut
# FIXME: ça marche pas
labelLayer(Selection_Communes_LR, txt = "NOM_COM", halo = TRUE, bg = "white", r = 0.05, cex = 0.6, pos = 3, font = 3, offset = 0)
Export de classif.
st_write(IRISDVFCAH, "data/IRISDVFCAH.gpkg", append = FALSE)
## Deleting layer `IRISDVFCAH' using driver `GPKG'
## Writing layer `IRISDVFCAH' to data source `data/IRISDVFCAH.gpkg' using driver `GPKG'
## Writing 67 features with 16 fields and geometry type Multi Polygon.
Zoom sur la commune de La Rochelle
# pour ne garder que la commmune de La Rochelle
cog_roch <- cog %>%
filter(nom == "La Rochelle")
# pour ne garder que les mutations correspondant
# aux appartements dans La Rochelle Agglomération (dans un premier temps)
dvf_appt_roch <- DVF %>%
filter(codecommune == 17300 & type == "Appartement")
# pour créer une liste d'années, dans l'ordre
dvf_appt_roch <- dvf_appt_roch[order(dvf_appt_roch$annee, decreasing = FALSE), ]
ListAnnees <- unique(dvf_appt_roch$annee)
# pour définir comme emprise pour le lissage, l'étendue
# de La Rochelle avec une zone tampon de 0,5 km autour afin
# d'englober les mutations limitrophes et réduire l'effet de bord
bbox_roch <- as.owin(st_buffer(st_geometry(cog_roch), 500))
# Création d'une palette de couleurs
cols <- c(
"#fff7f3",
"#fee1d3",
"#fcbea5",
"#fc9777",
"#fb7050",
"#f14432",
"#d42020",
"#ad1016",
"#67000d",
"#480000",
"#300000")
Série de cartes lissées montrant l’évolution des prix au m² des appartements calculés directement à partir des mutations de la Rochelle
# Paramétrage des marges pour insérer le titre général et les titres de chaque carte
par(
oma = c(3.5, 0, 3, 0) + 0.1,
mar = c(0, 0.5, 1.2, 0.5)
)
# pour découper la fenêtre en 3 lignes et 3 colonnes (9 années)
par(mfrow = c(3, 3))
# boucle pour produire et cartographier une surface lissée par année
for (i in ListAnnees){
print(i)
# Recuperation des jeux de données par année
IRISDVF_CA <- DVF %>%
filter(annee == i) %>%
drop_na(prix_m2)
print(paste("Number of mut =", length(IRISDVF_CA$prix_m2)))
print(paste("Number of mut =", mean(IRISDVF_CA$prix_m2)))
print(paste("Number of mut =", sum(IRISDVF_CA$prix_m2)))
IRISDVF_CA <- aggregate(
dplyr::select(IRISDVF_CA, prix_m2), by = IRISLR_CA, FUN = mean
)
IRISDVF_CA <- IRISDVF_CA[which(!is.na(IRISDVF_CA$prix_m2)), ]
# pour récupérer les coordonnées du centroïde des IRIS
pts <- st_coordinates(st_point_on_surface(st_geometry(IRISDVF_CA)))
# pour créer un objet ppp (format spatstat) et y intégrer dedans
# l'emprise et les valeurs à lisser (prix au m²)
IRISDVF_CA.ppp <- ppp(pts[,1], pts[,2], window = EmpriseLR, marks = IRISDVF_CA$prix_m2)
# pour calculer la surface lissée (rayon lissage : 250m et resolution
# spatiale de l'image : 100 m)
cartelissee <- Smooth(IRISDVF_CA.ppp, sigma = 250, weights = IRISDVF_CA.ppp$marks, eps = 100)
# Conversion de la surface lissée au format raster
cartelissee.raster <- rast(cartelissee)
crs(cartelissee.raster) <- st_crs(cog_roch)$srid # pour spécifier un SCR à l'objet raster
cartelissee.raster <- mask(cartelissee.raster, cog_roch) # découpage sur emprise La Rochelle
# Définition manuelle des seuils
bks <- c(min(values(cartelissee.raster), na.rm = TRUE), 1500, 2000, 2200, 2400, 2800, 3000, 3200, 3500, 4000, max(values(cartelissee.raster), na.rm = TRUE))
# reclassification de la surface lissée
cartelissee.reclass <- classify(cartelissee.raster, rcl = bks)
# vectorisation de la surface reclassée
cartelissee.vecteur <- st_as_sf(as.polygons(cartelissee.reclass))
# trier les valeurs
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$lyr.1, decreasing = FALSE), ]
# Tracer la carte
plot(st_geometry(cog_roch), border = "white", bg = "grey90")
typoLayer(
x = cartelissee.vecteur, var = "lyr.1", col = cols, lwd = 0.1,
border = NA, legend.pos = "n", add = TRUE
)
plot(
st_geometry(IRISLR),
border = adjustcolor("white", alpha.f = 0.5),
lty = 3, lwd = 0.01, add = TRUE
)
title(main = paste("", i, sep = ""))
}
## [1] 2014
## [1] "Number of mut = 2379"
## [1] "Number of mut = 2822.79118247799"
## [1] "Number of mut = 6715420.22311514"
## [1] 2015
## [1] "Number of mut = 2666"
## [1] "Number of mut = 2764.65239661806"
## [1] "Number of mut = 7370563.28938374"
## [1] 2016
## [1] "Number of mut = 3099"
## [1] "Number of mut = 2763.4410373335"
## [1] "Number of mut = 8563903.77469651"
## [1] 2017
## [1] "Number of mut = 3733"
## [1] "Number of mut = 2900.45902116142"
## [1] "Number of mut = 10827413.5259956"
## [1] 2018
## [1] "Number of mut = 3827"
## [1] "Number of mut = 3005.90376906455"
## [1] "Number of mut = 11503593.72421"
## [1] 2019
## [1] "Number of mut = 4093"
## [1] "Number of mut = 3216.62295294932"
## [1] "Number of mut = 13165637.7464215"
## [1] 2020
## [1] "Number of mut = 2199"
## [1] "Number of mut = 3457.52982466441"
## [1] "Number of mut = 7603108.08443704"
## [1] 2021
## [1] "Number of mut = 2409"
## [1] "Number of mut = 3702.98092072082"
## [1] "Number of mut = 8920481.03801646"
## [1] 2022
## [1] "Number of mut = 1251"
## [1] "Number of mut = 4000.58725714111"
## [1] "Number of mut = 5004734.65868353"
barscale(
lwd = 1.5, cex = 0.6, pos = "bottomright", style = "pretty", unit = "m"
)
cartography::north(pos = "topleft")
# Pour afficher le titre principal et la source
mtext(
"Prix moyen en Euros au m² des appartements à La Rochelle de 2014 à 2022",
cex = 1, side = 3, line = 1, adj = 0.5, outer = TRUE
)
mtext(
" Sources : IGN et DGFip - Rayon de lissage : 250 m\n et résolution : 100 m",
side = 1, line = 1, adj = 0, cex = 0.6, outer = TRUE
)
# Affichage de la légende (dans une nouvelle fenêtre graphique)
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot.new()
legendChoro(
pos = "bottomright",
title.txt = "",
breaks = bks,
nodata = FALSE,
values.rnd = -1,
col = cols,
border = "white",
horiz = TRUE
)